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Abstract. We report results of a study of the Newtonian dynamics of A' self-gravitating particles which start in a quasi- 
uniform spherical configuration, without initial velocities. These initial conditions would lead to a density singularity at the 
origin at a finite time when N ^ 'x>, but this singularity is regulated at any finite A' (by the associated density fluctuations). 
While previous studies have focussed on the behaviour as a function of N of the minimal size reached during the contracting 
phase, we examine in particular the size and energy of the virialized halo which results. We find the unexpected result that the 
structure decreases in size as A' increases, scaling in proportion to A'^'/^, a behaviour which is associated with an ejection of 
kinetic energy during violent relaxation which grows in propertion to A''/^. This latter scaling may be qualitatively understood, 
and if it represents the asymptotic behaviour in A' implies that this ejected energy is unbounded above. We discuss also tests 
we have performed which indicate that this ejection is a mean-field phenomenon (i.e. a result of collisionless dynamics). 
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INTRODUCTION 

We will discuss here numerical simulations, and some simple analytical results, of what appears at first sight perhaps 
a very simple initial condition which one might consider naturally in the study of the vmalization of self-gravitating 
systems: a large number of point particles are randomly distributed in a sphere, and evolved from a cold (i.e. zero 
velocity) start. The "simplicity" of the family of initial conditions is that they are characterized by a single parameter, 
the particle number A^ (as the system is open the units may be defined by the system size, its total mass and Newton's 
constant G). Further as N ^ °° the initial conditions tends to that considered as a starting point for analyses of non- 
linear structure formation in cosmology — the "cold spherical collapse" model, in which a perfectly uniform spherical 
overdensity embedded in an expanding universe is considered. A random sampling with particles of such a flat density 
profile is the simplest discrete realisation of this theoretical model one can envisage, and might thus be expected 
perhaps to be the subject of much study. In practice this initial condition has been but little studied — references to 
the few previous studies will be given below — because of the intrinsic difficulty in its numerical integration which is 
related to a property of the uniform limit: asN the evolution of the system leads to a singularity after a finite time, 
as all the mass arrives at the origin after a time t^ou ^ 1 / ^/Gpo where po is the initial mass density. While for any finite 
A^ system the singularity does not occur, the typical size of the region the system contracts to before "turn-around" 
decreases as A^ increases, and the typical particle velocities grow. This makes the numerical integration very costly, 
and hmits greatly the accessible particle A^ compared to other warm or less homogeneous initial conditions. 

One motivation for this study comes thus from the "uniform spherical collapse model": when calculating predictions 
for the masses and abundances of halos in the framework, one of the critical assumptions is that all the mass and energy 
in the initially collapsing region is ultimately virialized in the collapsed structure. The question arises as to whether 
this is generically true independently of the initial conditions. In the present case — which we have noted is, in a 
simple sense, the "closest" initial condition to the exactly uniform case — it turns out, interestingly, that this is not a 
good approximation: as we will describe, the violence of the collapse leads to an ejection of energy from the system 
(as kinetic energy of particles which escape with positive energy). Our numerical analysis, coupled to an analytical 
scaling argument, lead us to conclude that the ejected energy is in fact unbounded above as A^ increases, so that these 
collapses can be characterized as causing "explosions", with purely Newtonian gravitational physics. 

The initial motivation for our study of these initial conditions was, however, another one: this class of initial 



conditions poses issues about discreteness (i.e. N-dependent, non-coUisional) effects in an interesting way. The aim of 
A^-body simulations in cosmology (and indeed, most such simulations in astrophysics more generally) is to represent 
the colhsionless hmit, which means that the results of such simulations should be explicitly independent of (or 
depend very weakly on them at the A'^ considered). In practice this means that, if the dynamics is indeed colhsionless, 
it should be possible to define an appropriate extrapolation of which gives stable macroscopic results. In the present 
class of initial conditions it is clear that the simple extrapolation of N ^ °° described above does not give A^- 
independent results — as we have already indicated we find expUcit A'^ dependences of, for example, the final size 
and binding energy of the virialized structure. There are then two possibilities: either the physics of this ejection is not 
in fact colhsionless — which one would perhaps imagine might be the case as the particles do coherently contract into 
a very small region — or else the extrapolation in A^ must be performed in a different way. Through a careful study 
which we will summarize below we have shown that it the latter which is the correct explanation. The physical reason 
is that the A' dependence which appears in the final state when one extrapolates "naively" is in fact a result of the fact 
that A' controls the amplitude of the fluctuations about uniformity, which are the physically relevant quantities. The 
appropriate extrapolation to the coUisionless limit is that in which the particle number is increased, while keeping these 
fluctuations fixed. This extrapolation requires the introduction of an additional characteristic scale, which is the scale 
above which the fluctuations are kept (approximately) invariant. The existence of the colhsionless limit is, conversely, 
a result of the fact that the dynamics is only sensitive to fluctuations above some finite scale. The study of this problem 
for these initial conditions thus allow one to understand some of the subtleties which may be involved in practicaUy 
testing for collisionality effects in A^ body systems. This is in principle an important issue in cosmological simulations 
(see [1] for a detailed discussion). 

This proceedings essentially summarizes synthetically results reported in much greater detail in [2]. Because of 
space limitations we refer the reader to this paper for a fuller set of references on previous studies of virialization, 
giving here only the essential references, in particular to the few previous detailed studies of these specific initial 
conditions [3, 4, 5, 6]. 



COLD COLLAPSE THEORY 

We recall first the limit N ^ oo. The radial position r{t) of a test particle in an (idealized) exactly uniform spherical 
distribution of purely self-gravitating matter of initial density po and initially at rest (at time ? = 0) is simply given by 
the homologous rescahng 

r{t)=R{t)r{0) (1) 
where the scale factor R{t) may be written in the standard parametric form 

/?(^) = i(l+cos(^)) (2) 
K^) = ^(^+sin(^)), 

with Tscm ^ 32Gpo ' ^^^^ System collapses to a singularity, and physical quantities diverging. More 

specifically, taking ^ =n — e and expanding to leading order in e gives {t — Xgcm ) ~ , we have that R{t )^[t — Xscm\ 
and therefore the test particle velocities v{t), proportional also to the initial radius r(0), scale as v{t) ~ [f — T^em]"^^^ 

A' randomly placed particles in a spherical volume can be treated, up to some time and at sufficiently large scales, 
as a perturbed version of this uniform limit. Consider the approximation in which we treat the perturbations as if they 
evolve also in an infinite contracting system (i.e. neglect the effect of the boundaries on the evolution of the density 
perturbations). In the manner standard in cosmology (for the case of an expanding universe) one can then consider the 
fluid limit for the system and solve the appropriate equations perturbatively (see e.g. [7]). In the eulerian formalism 
this gives, at hnear order, a simple equation for 5(x), the density fluctuation (with respect to the mean density): 

5 + lH5-AnGpa5 = (3) 

where H{t) = R/R (dots denotes derivatives with respect to time) is the contraction ("Hubble") rate. These equations 
are derived in "comoving" coordinates x = r/R{t), where r are the physical vector positions. Note that /?x == r — Rx 
i.e., /?(?)x(?) is the "peculiar" velocity with respect to the "Hubble flow". 



It is straightforward to show, from Eq. (3), that, in the Umit/? ^ 1, 



(4) 

This is simply the usual decaying mode of the expanding EdS universe, which becomes the dominating growing mode 
in the contracting case. 

The singular behaviour of the spherical collapse is regulated by the fluctuations present at any finite in the initial 
conditions we study. A simple estimate of the scale factor at which one expects the spherical collapse model to 
break down completely may be obtained by assuming that this will occur when fluctuations at some scale (e.g. of order 
the size of the system) go non-linear. For Poisson distributed particles we have a mass variance proportional to A'^, and 
so the amplitude of the initial normalized density fluctuations is proportional to 1/ v^. Using the growth law given in 
Eq. (4), we can infer [3, 4]. 

Rmin^N-^^K (5) 

NUMERICAL RESULTS 

Details of our numerical simulations, performed using the publicly available and widely used GADGET2 code [8, 9], 
are given in [2]. Our study covers a range of between several hundred and several hundred thousand. We mention 
here just one important consideration: instead of the exact Newtonian potential, the code employs, for numerical 
reasons, a two-body potential which is exactly Newtonian above a finite "smoothing length" e, and regularized below 
this scale to give a force which is 1) attractive everywhere and 2) vanishes at zero separation. The (complicated) 
analytic expression for the smoothing function may be found in [9]. With this modified force the code does not integrate 
accurately trajectories in which particles have close encounters, which lead to very large accelerations and thus the 
necessity for very small time steps (which is numericaUy costly). However, on the (short) time scales we consider 
such trajectories should not play any significant role in modifying the macroscopic properties we are interested in. The 
results shown below correspond to a constant value of e in all simulations, a few times smaller than £ in the largest 
simulation. As we wiU detail further below when we discuss the coUisionless limit for our system, we have tested our 
results in particular for stability when e is extrapolated to smaller values, and we interpret them to be indicative, on 
the relevant time scales, of the e = 0, i.e., the exact Newtonian limit' . 

Qualitatively the evolution we observe in all our simulations is the same, and hke that well known in both 
astrophysics, and, more generally, in statistical physics for systems with long-range attractive interactions from sub- 
virial initial conditions of this type (i.e. with an initial virial ratio larger than -1) : the system first contracts and relaxes 
"violently" (i.e. on timescales of order the dynamical time scale Vm) to give a virialized, macroscopically stationary, 
state (see e.g. [10]). 



Minimal size: phenomenology 

The existing studies in the astrophysical literature of this class of initial conditions [3, 4] focus on how the singular 
coUapse of the uniform spherical collapse model is regulated at finite A'', and in particular on the scahng with N in 
numerical simulations of the minimal size reached by the system. Indeed in the study of [3] the "points" represent 
masses with extension (e.g. proto-stars) and the central question the authors wish to address is whether these masses 
survive or not the collapse of a cloud of which they are the constituents. This minimal radius Rmin may be defined in 
different ways, e.g., as the minimal value reached by the radius, measured from the center of mass, enclosing 90% of 
the mass. Alternatively it can be estimated as the radius inferred from the potential energy of the particles, the minimal 
radius corresponding to the maximal negative potential energy. The behavior of Rmin, determined by the first method, 
as a function of A^ is shown in Fig. 1. The fitted fine Rmin °= N~^^^ is the theoretical behavior predicted by the simple 
arguments given in the previous section. Agreement with this simple prediction has also been verified in both [3] and 
[4], the latter for an A'^ as large as 10^. 



' As a test we have also performed simulations using a code with a direct summation, and without any smoothing. For the range of N (up 
to a few thousand) for which we can run this code over the same physical time-scale, we find excellent agreement with the results obtained with 
GADGET2 with the smoothing we have adopted (see [2] for exact parameter values, as well as details of energy conservation etc.). 
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FIGURE 1. Behaviour of the minimal radius R„i„ attained, determined as described in text, as a function of N. The solid line is 
the best fit to the prediction of Eq. (5). 
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FIGURE 2. The behavior of /^(f), the fraction of the particles with positive energy, as a function of time for two different 
simulations. 



Mass ejection 

All particles start with a negative energy, but a finite fraction can in principle end up with a positive energy and 
escape from the system completely. While evidently the ejected mass is bounded above by the initial mass, the ejected 
energy is, in principle, unbounded above as the gravitational self energy of the bounded final mass is unbounded below. 

Shown in Fig. 2 is the fraction fP of the particles with positive energy as a function of time in two different 
simulations, while Fig. 3 shows the asymptotic value of in each simulation as a function of (i.e. the value attained 
on the "plateau" in each simulation after a few dynamical times, corresponding to particles which are definitively 
ejected on these time scales).While some previous works (see [2] for references) have noted the ejection of some small 
fraction of the mass in similar cases, the significance of the energy ejection as A' increases, and its dependence, 
has not previously been documented. Theoretical studies of the ejection of mass from a pulsating spherical system — 
which is quaUtatively similar to that described below for the ejection observed here — can be found in [1 1, 12]. 

Although /'' fluctuates in different realizations with a given particle number, it shows a very slow, but systematic, 
increase as a function of A^, varying from approximately 15% to almost 35% over the range of A' simulated. A 
reasonably good fit is given by 

r{N)^a + b\og{N), (6) 

where a = 0.048 and b = 0.022. Alternatively it can be fit quite well (in the same range) by a power law f « 0. IN^-^. 
Note that these fits cannot, evidently, be extrapolated to arbitrarily large A^ (as the mass ejected is bounded above), 
and thus our study does not actually definitively determine the asymptotic large A'^ behavior of this quantity despite the 





y = 0.048 + 0.022 ln(x) 
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FIGURE 3. Behavior of the fraction of ejected particles as a function of the number N of simulation particles in the system. The 
solid line is the phenomenological fit given by Eq. (6). 



large particle numbers simulated. As we will discuss briefly below, however, the mechanism we observe for this mass 
ejection leads us to expect that the value of fP should saturate when fP ~ 0.5. 



Energy ejection 

Let us now consider the energy carried away by these particles. Using total energy conservation, and the fact that 
both the final potential energy of the ejected particles and that associated with the interaction of the bound and ejected 
particles is negligible, we have 

Eo = W'' + KP + K" . (7) 

where £0 is the initial energy of the system, W" and K" are the potential and kinetic energy of the particles which are 
finally bound, and RP is the kinetic energy of the escaping particles. Further, since the bound particles in the QSS are 
virialized we have 

2K" + W" = 0. (8) 

Thus we have 

W" = -2K" = 2{Eo - RP) (9) 

Fig. 5 shows the ratio RP / fP, i.e., the kinetic energy per unit ejected mass, as a function of A^. Its behavior is fit 
very well by RP /fP N^l^ . Note that for the largest values of A'^ simulated RP is almost ten times the initial (potential) 
energy £0 of the system. It follows that we have the approximate behavior W" «: —N~^l^ (when we neglect the slow 
observed variation with of fP). 



Properties of virialized "halo" 

The dependence of the ejected energy on implies that the macroscopic properties of the virialized structure (which 
we refer to as a "halo" in the sense current in cosmology) also depend on A'. Studying the radial density profiles of the 
(approximately sphericaUy syimnetric) halos we find that they can always be fit well by the simple functional form 

n{r)= , "V,, ■ (10) 



This form of the profile agrees well with that found in previous studies for collapses from low initial virial ratio (see 
[2] for references). The A'^ dependence we find here is encoded in that of the two parameters n^ and ro, which we find 




FIGURE 4. Behaviour of the total kinetic energy for two simulations with different number of particles. 
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FIGURE 5. Observed behaviour of the ratio / as a function of particle number N. 



are well fit by ro °= N~^^^ and riQocN'^. In Fig. 6 we show the density profiles for various simulations with different A'^ 
where the axes have been rescaled using these behaviors. 

It is simple to show that these scalings with of no and ro are simply those which follow from those just given for 
f and W": using the ansatz Eq. (10) one has that the number A^" of bound particles is proportional to ne'e while the 
potential energy W" is proportional to m^ngrg where m is the mass of a particle. The best fit behaviors for ro and no 
thus correspond, since w 1 /N, to N" ^ N (i.e. a constant bound mass, and therefore a constant ejected fraction of 
the mass /'') and W" ~ N^^'^. More detailed fits to no and ro show consistency also with the very slow variation of 
/'' observed. In summary the dependence of the virialized structure manifests itself to a very good approximation 
simply in a scaling of its characteristic size in proportion to Rmin, the minimal radius attained in the collapse (which, 
as we have seen, is proportional to the initial inter-particle separation). 



FURTHER ANALYSIS AND DISCUSSION 
Mechanism of energy ejection 

We limit ourselves here to a very brief qualitative description of a detailed study we have performed (see [2]) of the 
evolution of the system during collapse and the mechanism which leads to the mass and energy ejection. It is simple 
to establish that the probability of ejection is closely correlated with particles' initial radial positions, with essentially 
particles initially in the outer shells being ejected. The reason why this is so can be understood as follows. Firstly, 




particles closer to the outer boundary systematically lag (in space and time) with respect to their uniform spherical 
collapse trajectories more than those closer to the center. This is an effect which arises from the fact that, when mass 
moves around due to fluctuations about uniformity, there is in a radial shell at the boundary no average inward flux of 
mass to compensate the average outward flux. The mean mass density thus seen by a particle in such a shell decreases, 
leading to a slowdown of its fall towards the origin. This "lag" with respect to particles in the inner shells propagates 
in from the boundaries with time, leading to a coherent relative lag of a significant fraction of the mass by the time of 
maximal compression. Secondly, these lagging particles are then ejected as they pick up energy, in a very short time 
around the collapse, as they pass through the time-dependent potential of the particles initially closer to the center, 
which have already collapsed and "turned around". This is iUustrated in Fig. 7, which shows, for a simulation with 
= 131072 particles, the temporal evolution of the components of the mass which are asymptoticaUy ejected or 
bound. More specifically the plot shows the evolution of Ve (and Vh) which is the average of the radial component 
of the velocity for the ejected (and bound) particles, and also e^ (and e,!,) which is the mean energy per ejected (and 
bound) particle (i.e. the average of the individual particle energies). The behaviors of Ve and v^, show clearly that the 
ejected particles are those which arrive on average late at the center of mass, with Vg reaching its minimum after the 
bound particles have started moving outward. Considering the energies we see that it is in this short time, in which the 
former particles pass through the latter, that they pick up the additional energy which leads to their ejection. Indeed the 
increase of e<. sets in just after the change in sign of v^,, i.e., when the bound component has (on average) just "turned 
around" and started moving outward again. The mechanism of the gain of energy leading to ejection is simply that the 
outer particles, arriving later on average, move through the time dependent decreasing mean field potential produced 
by the re-expanding inner mass. Assuming that the fraction of the lagging mass is independent of A^, an analysis of the 
scaling (see [2]) of the relevant characteristic velocity/time/length scales allows one to infer the observed scaling of the 
ejected energy with A^. Quantitatively we have not been able to explain, on the other hand, the observed dependence 
of the lagging mass, which should determine the A'^ dependence of /''. Given, however, that it is determined by a lag of 
the outer mass relative to that of the iimer mass, it seems clear that, as required, the mechanism observed will naturally 
lead it to saturate at a fixed fraction, of order one, as A'^ increases arbitrarily. 



Is dynamics coUisionless? 

The mechanism we have described for the mass and energy ejection is clearly of the "mean-field" type — the 
particles have been considered to move in the field produced by the bulk of the mass and no role has made ascribed, 
notably, to collisions with individual particles. As the system has contracted so much before re-expanding, particular 
care should evidently be taken in verifying that this is indeed the case, i.e., that the simulations do indeed represent well 
the coUisionless limit. As discussed briefly this is none other than a particular case of a question which can be posed 
about gravitational A^-body simulation, but one which shows the subtleties there may be in defining an appropriate 
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FIGURE 7. Radial velocity, and average energy per particle, as a function of time, of particles which are bound/ejected at the 
end of the simulation of 131072 particles. The energy of the particles has been arbitrarily rescaled. 



numerical extrapolation to provide a clear answer. More precisely the question is whether the system approximates 
well the evolution which would be obtained from a set of "collisionless Boltzmann equation" coupled to the Poisson 
equation, i.e. the Vlasov-Poisson (VP) equations. Is this the case? To determine whether it is we need to understand 
how we can test for its validity. As the VP limit is an appropriate N limit for the system, this means specifying 
precisely how this limit should be taken. We can then extrapolate our numerical simulations to larger N to test for the 
stability of results. 

Firstly it is clear that the appropriate extrapolation for the system we have studied is not the naive limit N ^ °°, i.e., 
in which we simply increase the number of Poisson distributed particles: we have explicitly identified macroscopic A' 
dependencies in fundamental quantities, so the results at any given A'^ do not approximate those at any other A'^, and 
indeed do not converge towards any A^-independent behavior. 

Formal proofs of the validity of the VP limit [13] for a self-gravitating system require, however, that the singularity 
in the gravitational force at zero separation be regulated when the limit N ^oois taken. This suggests we should take 
the limit N ^ °° while keeping fixed a smoothing scale, like the e we have introduced in our simulations. Doing so we 
would indeed expect to obtain a well defined A^-independent result, corresponding to the uniform spherical model with 
such a regularization of the force: the sphere will not collapse below a radius of order e, as the force is then weaker 
than the Newtonian force (and goes asymptotically to zero). One would then expect to obtain, for sufficiently large A^, 
a final state which is well defined and A'^ independent, but dependent on the scale of the regularization and indeed even 
on the details of its implementation. 

This limit is not the VP limit relevant here. We have indeed introduced a regularization of the force, but this has been 
done, as we have discussed, only for reasons of numerical convenience, and our criterion for our choice of e is that it be 
sufficiently small so that our numerical results are independent of it. Our results are thus, a priori, independent of the 
scale e (and of how the associated regularization is implemented). To illustrate this we show in Fig. 8 the evolution of 

(the fraction of particles with positive energy) as a function of time, in simulations from identical initial conditions 
with A^ = 32768 particles in which only the value of e has been varied, through the values indicated. Other quantities 
we have considered show equally good convergence as e decreases. Note that for the given simulation the mean inter- 
particle distance in our units is ^ = 0.016, so that the convergence of results is attained once e is significantly less 
than £. As we have seen, the minimal size reached by the collapsing system scales in proportion to £. We interpret the 
observed convergence as due to the fact that the evolution of the system is determined primarily by fluctuations on 
length scales between this scale and the size of the system. Once e is sufficiently small to resolve these length scales 
at all times, convergence is obtained. 

Given the essential role played by fluctuations to the mean density in determining the final state, it is clear that 
only an extrapolation of A^ which keeps the fluctuations in the initial conditions fixed can be expected to leave the 
macroscopic results invariant. Any change in N necessarily leads, however, to some change in the fluctuations. If, 
however, as indicated by the above results, there is a minimal length scale in the initial conditions for "relevant" 
fluctuations, we expect an extrapolation of A'^ which leaves fluctuations above this minimal scale unchanged to give 
stable results. 
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FIGURE 8. Evolution of the fraction of the mass with positive energy for simulations with N = 32768 for the different values 
indicated of the smoothing parameter e. 



Such an extrapolation for our initial conditions can be defined as follows. Starting from a given Poissonian initial 
condition of particles in a sphere, we create a configuration with A^' = nN particles by splitting each particle into 
n particles in a cube of side 2rs, centered on the original particle. The latter patlicles are distributed randomly in the 
cube, with the additional constraint that their center of mass is located at the center of the cube, i.e., the position 
of the center of mass is conserved by the "splitting". In this new point distribution, which has the same mean mass 
density as the original distribution, fluctuations on scales larger than are essentially unchanged compared to those in 
the original distribution, while fluctuations around and below this scale are modified (see [14] for a detailed study 
of how fluctuations are modified by such "cloud processes".). We have performed this experiment for a Poisson 
initial condition with N = 4096 particles, spUtting each particle into eight (n = 8) to obtain an initial condition with 
A^' = 32768 particles. Results are shown in Fig. 9 for the ejected mass as a function of time, for a range of values 
of the parameter r,, expressed in terms of £, the mean inter-particle separation (in the original distribution). While 
for rg = O.Si the curve of ejected particles is actually indistinguishable in the figure from the one for the original 
distribution, differences can be seen for the other values, greater discrepancy becoming evident as rj increases. This 
behavior is clearly consistent with the conjecture that the macroscopic evolution of the system depends only on initial 
fluctuations above some scale, and that this scale is of order the initial inter-particle separation £. And, as anticipated, 
this tianslates into an independence of the results when A^ is extrapolated in this way for an smaller than this scale. 

This prescription for the VP limit can be justified theoretically using a derivation of this Umit through a coarse- 
graining of the exact one particle distribution function over a window in phase space (see e.g. [15]). The VP equations 
are obtained for the coarse-grained phase space density when terms describing perturbations in velocity and force 
below the scale of the coarse-graining are neglected. A system is thus well described by this continuum VP limit if 
the effects of fluctuations below some sufficiently small scale play no role in the evolution. The definition of the limit 
thus requires exphcitly the existence of such a length scale, and the limit is approached in practice when the mean 
inter-particle distance becomes much smaller than this scale. With the kind of procedure given we have defined not 
only an extrapolation of N which gives stable results, but also a method of identifying this scale. 



Some conclusions 

The most surprising result of our study is that a simple initial configuration of particles interacting by Newtonian 
gravity only can liberate, in a time of order the dynamical time, an energy which is, apparently, unbounded above, 
growing approximately in proportion to N^^^. We have given a physical description of this ejection which accounts for 
this scaling, and performed careful numerical tests to check that the evolution is indeed coUisionless. This latter step 
required the definition of the appropriate numerical exttapolation, which in this case is quite non-ttivial. This provides 
an iUustration of the subtieties which may be involved in determining finite A^ corrections to the coUisionless limit, 
which we have argued elsewhere [1] are important notably in the context of cosmological simulations of dark matter. 
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FIGURE 9. Evolution of the fraction of particles with positive energy as a function of time, for the different indicated values of 
the parameter described in text. The "original" initial conditions has 4096 particles while the others have 32768 particles. The 
curve for = 0.8 is not visible because it is superimposed on that of the original distribution. 



In considering possible astrophysical relevance of these "classical gravitational explosions", a question which 
evidently arises is whether the energy ejection, which is associated essentially with the strong contraction during the 
collapse phase, is specific to initial conditions with spherical symmetry. In this respect we note that, although we expect 
smaller contraction factors for a generic initial condition, very large collapse factors may occur for non-symmetric 
initial conditions. This question has been considered at length in [4], which presents a detailed numerical study of the 
generalisation of the fluid SCM to axisymmetric and triaxial configurations. In the former singularities remain intact, 
and a relation Rmin °= N^^l^ is found empirically to replace the Rmin °^ N^^^^ behaviour of the spherical case. In the 
triaxial case the collapse factors are found to be typically finite, but they can be very large and no upper bound is 
placed on them. In forthcoming work we will address the ejection of mass and energy also for initial conditions with 
non-zero initial velocities. The divergence we have identified will be regulated again in this case, but this does not 
exclude that significant energy ejection may occur and be relevant, in particular, to understanding the properties of the 
"remnant" virialized structure. 
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